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far from equilibrium and, in particular, we find the exact probability distribution of 
the energy current in both the high- and low-temperature phase. 
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1. Introduction 

The infinite range Ising model, in which each individual spin interacts with the remain- 
ing — 1 ones has served as a useful testbench for many ideas in various subfields of 
statistical physics, ranging from critical dynamics to spin glasses. The reason is twofold: 
it is relatively easy to come up with exact yet nontrivial results for this system (in the 
large system limit at least), while at the same time it stands for a genuine interacting 
system which further possesses a phase transition from a disordered into an ordered 
state as the temperature is lowered. In the present work we will follow Ruijgrok and 
Tjon jl] who, by endowing the system with spin-flip dynamics, provided the first exam- 
ple of an exactly solvable critical dynamics problem back in 1973. In the meantime the 
subjects of interest have drifted towards other issues, but the technical motivations for 
using the infinite range Ising model remain. Our physical motivations have their roots 
in the study and identification of the generic properties of nonequihbrium steady-states 



Among the very few exact statements that can be made on NESS, the recent fluc- 
tuation -or Gallavotti-Cohen- theorem |2j plays a prominent role. This theorem is a 
symmetry property of the entropy current distribution function. As such it provides, 
at least formally, a prescription for obtaining an infinite set of Green-Kubo like rela- 
tions connecting current fluctuations to the system's response upon external forcing. 
While establishing a Gallavotti-Cohen theorem for Markov processes can be achieved 
in general terms 01^, explicit computations of current distribution functions are rare. 
For systems close to equilibrium such as boundary-driven lattice gases, Bodineau and 
Derrida [HI, and Bertini et al. \Q\ have provided a general method for computing cur- 
rent distribution functions from the sole knowledge of the Onsager response coefficients. 
However for the generic case of systems maintained in a steady-state far from equilib- 
rium, no general principle has hitherto been unraveled. To the best of our knowledge a 
single exact calculation exists for the totally asymmetric exclusion process, a degenerate 
case for which the Gallavotti-Cohen theorem (in the version of Lebowitz and Spohn 
does not hold. For this reason we have sought to exhibit a system driven far from 
equilibrium for which -albeit mean-field- such a current distribution function is accessi- 
ble: an Ising model in which spins are connected to heat baths at different temperatures. 

There are several ways of driving the infinite range Ising model into a NESS. One 
that is inspired from a series of recent works [3 El El consists in coupling spins to 
independent heat baths, thus creating a macroscopic energy current by means of a bulk 
drive. In the particular version we have coined, N Ising spins (Tj have a ferromagnetic 
interaction energy 



(NESS). 




(1) 
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and the dynamics of N/2 spins is generated by a heat bath of inverse temperature /3i 
while the dynamics of the remaining half of the spins is driven by another heat bath 
at a different inverse temperature 132- This state of affairs, namely the existence of two 
heat baths at unequal temperatures, leads to a steady energy current flowing through 
the spins from the warmer bath towards the colder one. 

The physical results we have obtained are concerned with the steady-state measure 
and the energy current distribution function. We have been able to provide an exact 
solution to the Fokker-Planck equation governing the probability distribution of mag- 
netization fluctuations, thus leading to the first example of an A^-body nonequilibrium 
system for which probability flow lines in phase space can actually be visualized. Our 
second achievement is to have provided the large deviation function of the energy cur- 
rent. The methods we have resorted to rest on the various formulations of the master 
equation governing the microscopic dynamics. On one hand, by means of a Van Kam- 
pen PI] expansion of the magnetizations around their mean values, we have obtained 
solvable Fokker-Planck equation describing steady-state fluctuations. Technically, this 
amounts to finding an eigenfunction of the master equation evolution operator. On the 
other hand, by means of a mapping of the master equation onto a quantum hamiltonian, 
we have determined the energy current distribution function, which, in technical terms, 
has proved to be an eigenvalue problem. 

Our first task will be to provide an accurate description of the stationary state 
distribution (Sec. 3), which we will present after having precisely characterized our 
model and its steady-state properties (Sec. 2). Sec. 4 will be devoted to a complete 
study of the energy current distribution function, in the light of the Gallavotti-Cohen 
theorem. Our concluding remarks will be followed by two appendices generalizing the 
aforementioned results to an infinite-range Ising model in contact with more than two 
heat baths. 

2. A two temperature Ising model: phase diagram and steady-state 
properties 

2.1. Microscopic dynamics 

The energy of a configuration of Ising spins <t = {(Xj} is given by 



We now divide the N spins into two groups with labels 1 and 2 of N/2 spins each. A 
spin cTj from set 1 flips with a rate 




(2) 



1=1 




(3) 
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Spins from group 1 try to equilibrate at inverse temperature /3i with respect to Ti. 
Similarly, a spin from group 2 flips according to 

Vj G 2, w,{a, ^ -a,) = e"/^^-^*^/^ (4) 

This is the infinite-range counterpart to the one-dimensional systems considered by Racz 
and Zia [7j, and Schmittmann and Schmuser [H]. Denoting by 

we see that when the temperatures are equal, Pi = P2 = P , or e = 0, the system reaches 
equilibrium at temperature p. This is because the rates fj3l4|) then satisfy detailed bal- 
ance with respect to the Gibbs distribution Z~^e~^'^. Though the precise expressions 
of the rates we have chosen differ from the original Glauber rates, they possess the 
same quahtative properties with some advantages in the large-system limit discovered 
by Ruijgrok and Tjon [T]. 



2.2. Phase diagram 

Introducing the mean magnetizations 



= J^C^^^j)^ "^2 = ^(Xl^i)' m = 1711+7712 (6) 

we may find the following evolution equations for the averages: 

— = — 2mi cosh PiTTi + sinh Pim, — — ^ = — 2m2 cosh P2m + sinh P2'm(7) 
at at 

from which one deduces that in the steady-state (provided it exists) 

m = - (tanh Pim + tanh P2'iti) (8) 

Interestingly, although the transition rates ()3l4j) are different from the standard Glauber 
rates, they lead to the same steady-state average magnetization. From (jH)) we deduce 
that in the steady-state the system undergoes a second order phase transition from a 
high-temperature disordered state at /3 < 1 in which mi = m2 = m = to a low- 
temperature ordered (doubly degenerate) state at /? > 1 with nonzero magnetizations. 
In the /? — >■ 1^ limit at e fixed one finds 

V3 r-r — - l + e VS r- — - 1-e VS 



m 



V/3-l> rm - ±— ^/ ^ Jp-1, 7712 ^ ±— JL—_ ^p-l{9) 



Vl + 3e2 ^ 2 VI + 3£2 ^ 2 VI + 3^2 

According to the magnitude of the nonequilibrium drive e it may be seen that the 
ordered state may be either ferromagnetic {\e\ < 1) or antiferromagnetic {\e\ > 1, if one 
allows for negative temperatures, as we shall discuss in our conclusion in Sec. E))- 
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2.3. Entropy and energy currents 

Following the prescription of Lebowitz and Spohn we may define a time integrated 
instantaneous entropy current by 

where o-(°) = o-(O), = cr(t) is the sequence of states occupied by the system over 
the time interval [0,t] (this is the history of the system between and t). The rates 
W{cr ^ a') of hopping from configuration cr to configuration cr' between t and t + dt 
are easily deduced from ()3I4|) . Inserting the explicit expressions for the W{a cr') 
leads to 

Qsit) = -/5 merit)] - n[a{Q)]) + sQ{t) (11) 
where we identify Q as the integrated energy current: 

n=0 

where aj^ is the spin being flipped at time n and M„ is the total magnetization at that 
moment. The sign + (resp. — ) corresponds to flipping a spin from group 1 (resp. 2). 
Note that TC[a(t)] — H[cr{0)] being bounded over time, Qs(t) and eQ(t) have the same 
large deviation functions. It is clear that on average, 

J, = Ml) = -l(5^a^.(M-a,)e-/^-^^^/^-5^a,(M-a,)e-^-^^^/^)(13) 

While interpreting Qs{t) as an integrated entropy current requires an elaborate 
reasoning the physical meaning of as an energy current is much more intuitive. 
Indeed, the total energy of the system is constant on average in the steady-state: 

^ = = -(Ji + J2) (14) 

where Jq, is the energy flux due to spin-flips caused by heat-bath a (for instance Ji is 
the first term on the rhs of (fT^ ). The quantity = Ji — J2 is therefore a measure of 
the energy flowing from group 1 towards group 2. Hence the related entropy current J5 
must read 

Js = PlJl+P2J2 = eJe (15) 

This interpretation of eJe as an entropy current has been discussed, on the grounds of 
phenomenological thermodynamics, by Racz and Zia [7]. 

As described in appendix B, there is no immediate link between the entropy current 
and an energy current for a system in contact with more that two heat baths. 



3. Stationary state distribution 

3.1. Van Kampen expansion and Fokker-Planck equation 

In this section we derive a Fokker-Planck equation governing the probabihlty P{xi,X2, t) 
of observing the following fluctuations of the spin magnetizations: 

Xa = — 1= , a = 1,2 (16) 



This is the Van Kampen |llj expansion of the master equation around the mean 
magnetizations nia- The \fN rescaling is precisely designed for the Xa to have order 1 
fluctuations. We find that P(xi,X2,t) satisfies the following Fokker-Planck equation: 

dtP = -d.,,Ji-d^,J2 (17) 

where the probability current is given by 

Ja = faiXuX2)P- D^d^^P (18) 

The two-dimensional force (/i, /2) does not derive from a potential unless both heat 
baths are at the same temperature. The general expression of the force components is 

fi{xi,X2) = ((/3i — 2)xi + l3iX2) cosh Pim — 2/5i(xi + X2)mi sinh Pim 
/2(xi, X2) = {{l32 - 2)x2 + l32Xi) cosh/52'm - 2l32{xi + X2)m2 sinh /52m (19) 

The diffusion constants are given by 



Da = cosh /3am — 2ma sinh /^^m = a/1 — 4m^ (20) 

In the high temperature phase, using that P1/2 = P ± e, this may easily be cast in the 
following form 

fi{xi,X2) = -d^JJe + eX2, f2{xi, X2) = -(9x2^4 - £Xi (21) 
where the potential energy has the expression 

TT ^ ^ n nAxl+X2f , iXl-X2f X\-X\ 

[/,(xi,X2) = (1 -/3) ^ + ^ e — ^ (22) 

In the high temperature phase, to which the ensuing analysis will be confined for 
simplicity (a general solution is provided in appendix A), where m\ = m2 = 0, we 
thus have to solve 

d^P = = _d^^j^_d^^j^ (23) 
where the probability current reduces to 

Ji = {{P-2)xi+pX2+e{xi+X2))P~d,,P, J2 = {{P-2)x2+pxi-e{xi+X2))P-d,,P{24) 
When Pi = P2 = P, that is in equilibrium, the distribution reads 

{Xi+X2)'^ (Xi-X2)^ 



P{xi, X2) ~ exp 



-a-py 



(25) 
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We may find the exact solution to the Fokker-Planck equation by having the intuition, 
following [TT], that the effective potential U^g defined by 

P{xi,X2) = Z-'expi-U^g) (26) 

will be quadratic in terms of Xi and X2- This is suggested by the force being linear in 
Xi and X2. And indeed this naive assumption leads to the effective potential U^ff given 
by 

U^n{xi,X2) = + X2f + {x^[l - J/2] - X2[l + J/2]f (27) 

2 4 + 

where we have introduced the constant J = ^^g- It is then an easy task to compute the 
mean energy current J^: 

J, = (^(-2(T|(Af^ - (T|)/Ar)e-^i"I^^/^ - J](-2(t|(AF - a;)/Ar)e-^2"I^'/^) 
= 2e{{xi+X2f)-2{{xl-xl)) 

In terms of the magnetization fluctuations the solution to the Fokker-Planck equation 
reads 



P{xi,X2) ~ exp 



^ ^(Xl + X2f - {x,[l - J/2] -X2[l + J/2]f 



(29) 



2 ^ ' 4 + 

A few comments are in order. In spite of the phase space being only two-dimensional, 
this is just enough to allow for inhomogeneous currents to flow (contrary to a one- 
dimensional phase space). While the total magnetization has global fluctuations equal 
to those of a system in equilibrium at /?, it may be seen that the magnetization difference 
between the two spin groups is increased with respect to its equilibrium counterpart in 
the presence of a current: 

((Xi - X2Y)j^o - {{Xi - X2)^)oq,j=o = 

This provides an example of a nonequilibrium drive giving rise to an increase of 
fluctuations, rather than to a decrease (as is usually noted, e.g. in driven lattice 
gases [121 Uni and spin chains |I1]). We now turn to an analysis of the probability 
flow lines. 



3.2. Flow lines 

In equilibrium, by definition, there is no probability current, while in a NESS there are 
steady (probabihty) currents. The flow lines, namely the set of points (xi, X2) such that 
Ji(xi,X2)dx2 — 72(2;!, 3^2) dxi = 0, in phase space turn out to be ellipses, as shown in 
Fig (0). For J 7^ the flow lines are ellipses of equation 

(1 -P)(l + ^] (xi + X2)^ + (xi[l - J/2] - X2[l + J/2]f = C (31) 
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Figure 1. The {xi,X2) axes denote the deviations from the group 1 and 2 
magnetizations. Flow lines of the probability current J' are represented for different 
values of the constant C in The vertical axe is the probability P{xi,X2), 

illustrating that the flow lines coincide with the isoprobability contours. On this figure, 
P = 0.8 and £ = 0.7. 



and they coincide with the isoprobabihty contours, a generic property of hnear force 
driven systems. For an equihbrium system the flow hnes can be seen to collapse onto 
a single point. Interestingly, the shape of the flow lines can be used to infer properties 
of the steady-state distribution: the departure from the ellipses will indicate deviations 
from linear forces in the Fokker-Planck equation. 



3. 3. Master equation and effective free energy 



We characterize a state of our system by the "local" magnetizations = X] 
and we denote (in this paragraph only) by rria = Ma/N. The steady-state solution to 
the master equation being denoted by P^^{Mi, M2), we define the effective free energy 
/(mi, ma) by 

lnP,t(A^mi, A^ms) 



/(mi, ma) 



lim 



(32) 



We split / into an entropic contribution s(mi,m2) and an effective energy e(mi,m2). 
Whether in equilibrium or in a NESS, the entropic part is defined by the combinatoric 
factor for the number of configurations with Mi and M2 : 



s(mi,m2) 



N 



In 



2 

jV+2Afi 
4 



/ K 
2 



yN+2M2 j 

1 + 2ma 1 + 2mc 
m ■ 



+ 



1 

m — 



2mn 



(33) 
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This fully defines e = f + s. In equilibrium we have that e^^{mi,m2) = — ^(mi + m2)^, 
and we wish to find how the nonequilibrium drive modifies this result, namely what kind 
of effective interactions between the two groups of spins it generates. Since ~ e^*^*"*^^ is 
a stationary solution to the master equation governed by the rates ()3l4p . we find that 
e(mi, 7712) is a solution to 

Q = (^l — 2mi) (^e"^i('"i+™2)-29ie _ g+/3i(mi+m2)j 
+ (1 + 2mi) (^Q+P'^^:^i+^-2)+'^9ie _ g-/3i(m.i+m2)j 

+ (1 — 2m2) (^e"^^*-™^'*'™^-'"^^^'^ — e"'"^2{mi+m2)j 

+ (1 + 2m2) (e+^2(mi+m2)+2a2e _ g-/32(mi+m2)j ^34) 

To first order in e one may verify that 

e(mi, 7712) = — l3m^ ""^("^i ~ m2)h{m) + 0{e^), m = mi + m2 (35) 

where the function h is the solution to the following first order ordinary differential 
equation: 

(m - tanh /3m) + /i(m) - (2 - /3)m = 0, /i(0) = (36) 
dm 

For instance, as m — » 0, 

/i(m) = m - — m^ + Oini^) (37) 

^ ^ 3(1 + 3(1-/3)) ^ ' ^ ' 

thus recovering the leading term of the high temperature Van Kampen expansion. At 
this stage we have completed our description of the steady-state properties. Note that 
the structure of ()35j) has flavors of the much more complex one found by Derrida et 
al. ^3] in the framework of the asymmetric exclusion process for the effective free energy 
of a given density profile. 



4. Energy current distribution 

This section is devoted to determining the large deviation function of the time integrated 
energy current. 

4.1. Modified master equation 

Let Q{t) be the fluctuating energy current integrated over the time interval [0,t] as 
defined in (fT2j) . We are interested in p{Q,t), the probability that Q{t) = Q at time t or 
in its generating function p{X,t) = {e~^^). It is possible to write a master equation for 

P{Mi,M2,Q,t) = Prob{^(7j- = Ml and ^a^- = M2 and Q{t) = Q}(38) 

Inserting the explicit expressions of the transition rates we arrive at 

N/2 + Ml + 2 M+2 , M + 1 N/2 - Mi , . m , 

dtP = \ e-^^— P(Mi + 2, M2, Q - 2^^, t) - ' ^ ^ +^^i^P(Mi, M2, Q, t) 
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+ iVi^ii±^,...M^P,M, - 2, M,,Q + 2^1^, - ^^^e-.#P(M., M,. Q, 

iV/2 - M2 + 2 , . A/-2 , ^ M - 1 , + M2 « M , ^ , 

+ ^ j;^ e+^^ — P(Mi, M2 - 2, g + 2^^, t) - ' ^ \ -^^^P{M^, M2, Q, t) 

(39) 

Going to the generating function 

P(Mi, M2, A, t) = ^ e-^^P(Mi, M2, Q, t) (40) 
Q 

and setting |^(A, t)) = Emi,M2 ^2, A, t) |Mi, M2), we may rewrite Eq. (jSHl) in the 

following form 

"^^-mm^i)) (41) 

where the operator ll{X) reads 

^(A) = 5^ (1 - aJe+^^'^.^(^^^-'^.^)/^) 6-^^^^!^^/^+ J] (1 - aJe-^^'^K^^^-'^l)/^) e-'^^'^.^^^/^(42) 

The asymptotic behavior of 

p(A,t) = ^ P(Mi,M2,A,t) = (p|e-^(^)*|^(0)), (p| = J] (Mi,M2| = projection state(43) 

will be governed by the largest eigenvalue /i(A) of —H{\) in the sense that 
lim^^oo jliip(A,t) = /i(A), which we now set out to determine. Before embarking into 
technicalities it is convenient, but by no means compulsory, to perform a similitude 
transformation on H{X) 

= ^ ^g-zSKT^MVAf _ ^a.g+(2A-e)a|MVA^-(2A+/3)/Af j 

+ ^ j^g-/32fr|MViV _ ^Xg-(2A-£)(r|MV7V+{2A-/3)/iV^ ^^^^^ 
i62 

The transformation (j44p does not have the same effect as that conducted by Ruijgrok 
and Tjon P -it does not make the resulting operator Hermitian- but it serves the 
same practical purpose: calculations are performed in a more convenient way where the 
system symmetries (upon exchanging the roles of 1 and 2) are made obvious. In terms 
of its symmetrized counterpart Hs{X), we have that 

HsWV = H,is - \*) (45) 



An important consequence of symmetry ()45|1 is that, for A real, both Hs{X) and Hs{6 — X) 
have the same spectrum, hence 

/i(A) = fx{e - A) (46) 
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This is the Gallavotti-Cohen theorem. A direct consequence for the energy current large 
deviation function 7c{q) = jlnp{Q = qt,t) is that 

7r(g) - 7i{-q) = eg (47) 

where we have used that 7r(g) = ma.xx{^{\) + \q}. Another useful consequence of 

is that for A G I + iR, is Hermitian, which will justify diagonalization in that 

region of the A complex plane. 



4-2. Mapping to a free boson problem 

We introduce, following Ruijgrok and Tjon pp, bosonic operators a^al^ (a = 1,2) to 
describe magnetizations 1 and 2 in the vicinity of the paramagnetic state: 

M: = N/2-2aia., My = -r^Nj2{al-a^), = v/iV72(al + a J (48) 

The relations (jlH|l hold provided we are interested in states such that the number 
operator aj,a„ remain of order unity, that is much smaller than \/N. In terms of these 
operators we find that 



i\ «! a\ 02 ) r(A) 



ai 



A)) (49) 



V a2 y 



with 



r 



V 



Z-2A 
Z + 2 
Z 

Z-{2\- 



Z + 2 
Z + 2{\-e) 
Z + 2\-e 
Z 



Z 

Z + 2\-e 
Z + 2\ 
Z + 2 



Z - (2A - 
Z 
Z + 2 
Z-2{X- 



(50) 



0/ 



where Z = — 1/5(2 — /5) + 2A(£: — A). For normalization purposes JH] it is necessary to 
define an auxiliary matrix F built from F by equating to zero in the latter all matrix 
elements connecting two creation or two annihilation operators. It is then a simple 
matter [12] to determine not only the ground-state but also the spectrum of -ffs(A). To 
do so we introduce the matrix Q defined as 

/ iu \ 




n 



\ 

















-ILO 



We must now evaluate the quantity 

det(f 



MA) = - /-In 

' 2 J 271 



^]) 



(51) 



det(F + ^]) 



-(/3^ + 4A(.-A)) 



(52) 



Given that 



det(F + fi) = (cu^ + ujDiu"^ + ul) 



(53) 
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with 



uj±{X) = V(2-/?)2 + 4A(£-A)± v//32 + 4A(£ - A) 
and that similarly 

det(f + fi) = {uj^ + cD+)(cu^ + c^!) 

with 



(54) 
(55) 



u^iX) = ^ ((2 - + 4A(e - A) ± ^[(2 - /5)2 + 4A(£ - A)] [/^^ + 4A(e - A)] 



we arrive at the following result: 



/i(A) = ^iCu+ + Cu_-u;^- -\{fi^ + 4A(e - A)) 



(57) 



which simplifies into 



^(A) = 2 - /5 - v/(2 - /9)2 + 4A(£ - A) (58) 

By the same token we obtain the spectrum of H{\), whose eigenvalues are given by 

Sp(i^(A)) = {Lu4\)i + u;4X)e - fxi\)}e,e'eN (59) 

Specializing to A = we obtain as a side result the spectrum of the master equation 
evolution operator, whose slowest relaxation time is given by uj^^{0) = (2(1 — 
This again matches the results of PJ. 

Given that 7r(g) and /i(A) are the Legendre transforms of each other we arrive at 
the explicit form of the current large deviation function 7r(g) 

n{q) = |g + 2 - /3 - + (60) 

and it has the graph shown in Fig.|2l 



7r{q) 



1.25 1.5 1 




Figure 2. This is tlie plot of tiie energy flux large deviation function 7r(g) (vertical 
axis) as a function of q (horizontal axis) at /3 = 0.5 and e = 0.4. 
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4-3. Energy current in the low temperature phase 



Similar methods allowed us to express the generating function of the cumulants of Q{t) 
in the low temperature phase, at /3 > 1. The final result is 



C1+C2-- — + — 
2 V ci C2 




Cl + C2 



1 M ^/52 

2 I Cl C2 



+ 



4 



C1C2 



■A(£- A)(61) 



where Cq, = cosh /3a (mi +7712) = 1/a/1 — 4m^ and where is the stationary solution of 
(|7j), and this is a function of /3 and e. However the current is not the order parameter of 
the phase transition, therefore nothing dramatic is expected to occur for /x(A) at /3 = 1. 
An important difference with the high temperature result must be emphasized: in the 
low temperature ordered regime the current is a nonlinear function of e, as plotted in 
Fig. 121 The similar mathematical structure of fi{X) in the high and low temperature 



Current J 
0.3^ 

0.25 

0.2 
. 15 

0.1 
. 05 



0.5 



1 . 5 



2 £ 



Figure 3. Plot of the average energy current J as a function of e € [0, /?] at /3 = 2 in 
the ordered phase. 

phases seems to be generically related to Langevin equations with linear forces [T71. 

The energy current at fixed drive e = 0.5 as a function of /? G [0.5, 2] represented in 
Fig. m shows that, from the disordered to the ordered phase, the current remains finite 
and continuous, though it develops a cusp at the critical point P = 1. 

4-4- Green-Kubo relations 

Exploiting the explicit formula (fSH|) for /i(A) we find, after differentiation with respect 
to A once and twice, that 

t 2-/3' t 2-/5^(2-/3)3^ ^ ' 
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Figure 4. Plot of the average energy current J as a function of /3 € [0.5, 2] at e = 0.5. 



Note that defining the diffusion coefficient D{f3) as the response to an external drive 
and (t(/3) as the variance of the current fluctuations we find 

_L^..(,).M.,^^ (a3) 

n=P2=p ^ — p r 



With these expressions one may verify an integral formulation of the Green-Kubo 
relation 

= e (64) 

Nevertheless the sole knowledge of -D(/3) and cr{[3) does not give access to the full 
distribution /i(A), as opposed to the cases studied by Bodineau and Derrida j3] by 
means of an additivity principle or by Bertini et al. [iHl E] who resorted to fluctuating 
hydrodynamics pT?| E] . In order the latter approaches to hold, the typical current must 
scale to zero with the system size at fixed (intensive) external field. This the second 
example, aside from the extensively studied Asymmetric Exclusion Process ^H], of an 
interacting system, albeit mean-field, in which the entropy (or energy) current can be 
computed exactly, with the additional property in our case that the Gallavotti-Cohen 
theorem is fulfilled, hence generalized Green-Kubo relations as well. 



5. Final comments 



We have been able to present explicit and exact results for the steady-state of a system 
made of interacting spins driven far from equilibrium by heat baths at different tempera- 
tures. The system described exhibits a ferromagnetic-to-paramagnetic phase transition. 
The simplicity of some of our results, like Gaussian fluctuations for the magnetizations, 
are admittedly an artifact of our infinite-range, mean-field, model. Nevertheless, due to 
easier mathematics, we have been able to precisely describe the probability flow lines, 
ellipses in a two-dimensional phase space. Other concepts arising within the framework 
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of dynamical systems theory, like that of topological pressure, once transposed to our 
model, might equally lend themselves to analytical approaches. 

Our other result of interest concerns the computation of an energy current dis- 
tribution for a system far from equilibrium, that cannot be described by fluctuating 
hydrodynamics, although it falls within the scope of the Gallavotti-Cohen theorem for 
Markov processes |3]. To our knowledge, this is the first one of this sort. It is only a 
first step towards the desirable, but remote, goal of characterizing stationary systems 
driven very far from equilibrium. 

Among the prospects uncovered by the present work, we mention the extension of 
the urn model of Bena et al. jSHj analyzed in the light of Greenblatt and Lebowitz' com- 
ment |2I]. In spite of being genuinely nonequilibrium, our model will most probably not 
display any surprises as far as the Yang-Lee zeros of the partition function are concerned, 
simply because the phase transition that takes place at /3 = 1 is akin to its equilibrium 
counterpart (and belongs to the same universality class). However the urn model may be 
very well defined for negative temperatures. Preliminary studies indicate that such rates 
open the door to limit cycles and chaotic behavior that we shall explore in future studies. 
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Academy of Sciences (Grant No. OTKA T043734). It is a pleasure for the authors to 
thank Henk Hilhorst, Cecile Appert and Bernard Derrida for their comments during the 
preparation of this work. 

Appendix A: Fokker— Planck equation for n heat baths 

In the disordered phase, the Ising system is split in n parts of magnetization fluctuation 
(1 < a < n), each of them equilibrated with a bath at inverse temperature 
/3q = + where (3 is the average of the /?a's. The Fokker-Planck stationary equation 
for the whole system can be written daJa = (we will use summation convention over 
repeated indices). The jT^'s are the components of the probability current, of expression : 



This is a special case of a Fokker-Planck equation with linear forces fa-yX^. As / is 
a definite negative matrix, we know from general considerations pT that its stationary 
state is Gaussian of covariance matrix : 



This matrix however proves difficult to be explicitated in the general case. Here we will 
restrict our analysis to the case where / can be diagonalized, as in our system, and use 
an indirect method to compute the result of (j(i7|l . First note that if / is symmetric. 



Jc 
fa 




(65) 

(66) 




(67) 
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setting U{xi, . . . ,x„) = —^x'^fx, we have fa-yX^ = —dJJ- In other words, the forces 
apphed on the system are conservative and deriving from the potential U . In that case, 
we immediately find the stationary solution of the Fokker-Planck equation by requiring 
Ja = ^ '■ this is, as expected, the Gibbs-Boltzman distribution P{x) ~ which 
corresponds to an equilibrium situation, and the defining absence of a probability cur- 
rent. 



For / diagonalizable but not symmetric, no solution can be found imposing Ja = 0. 
Our aim is to parallel the resolution of the equilibrium case, by transforming the Fokker- 
Planck equation into some new equation daJla = whose solution can be found by 
requiring = 0. We first perform the change of variable y = bx, where 6 is a matrix 
such that bfb^^ = Diag(Ai, . . . , A„). The Fokker-Planck equation takes the form : 

- baa'b^a'dadyP + <9a(A„?/„P) = (68) 

where now da = d/dt/a- The force term is symmetric, and can be formally written as 
deriving from some potential. There are many ways to split the first term of so as 
to write it as a divergence. Let A^-y be arbitrary constants. Writing, for a < 7: 

d^P) + d^{{l - Aa^)baa'b^a'daP) , (69) 

the FP equation now takes the announced form daJ^'a = 0, and requiring J''^ = means 
(defining ba ■ b^ = baa'bya') '■ 

( bl (2Aa,K ■ b^ \ 



hi 



VP = Diag(Ai,...,A„)y (70) 



\^ A^a)ba-b^^ bl j 



L 



This equation has a solution if, and only if A is symmetric, which occurs when: 

A., = (71) 
In terms of the original x variables, we finally get the following expression : 

P{x) = exp [-Ix'^b'^M-^bx] with (Mj^)„^ = f ^° " (72) 
\ 2 / ■ Afj + A^ 

in which we can read the result of (j67p in the general case. When applying this result 
to the n baths Ising problem, we find : 

1t,.-i\ ........ 1 / 1-/5 , {l+ea){l+e,) 



P{x) = exp --x' M~'x with : Af™ = -5™ H — ^ -7:^ — ^ + 



2 ; n ' n2 V 2-/3 ' (l-/3)(2-/3) 

The method above can be generalized to FP equations whose diffusion term is Tajdad^P 
(as, for instance, in the low-temperature phase of our Ising model). The final result takes 
the form (f7^ with {My)a^ = 2Ta''y'baa'b^^' / {\a + A^). 
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Appendix B: entropy current distribution function for n heat baths 

As in appendix A, we consider an Ising system split in n parts equilibrated with baths at 
inverse temperature (3a with 1 < a < n. The Lebowitz and Spohn 4j integrated entropy 
current (|Tn|l cannot be decomposed as in (|TT| . We will thus study the distribution 
of Qsif)- As in section HI the corresponding function //(A) is given by the largest 
eigenvalue of some operator —H{X), which reads : 

n 

H{X) = J2Y1 (e""^'^'^*'^^'^ - aje(2^-^)^"'^l*^//^) e"^'^/^ (74) 

where we defined the operator = ^i^jC^l verifying [MJ,ctJ] = 0. As expected on 
general grounds E] , this operator possesses the Gallavoti-Cohen symmetry : for A 
real, (^f(A))^ = — A). There is no need here to perform a symmetrization analogous 
to because we directly focus on the toto/ entropy current Qsif)- 

For the sake of simplicity, we will now confine our analysis to the disordered phase 
(/? < 1). In that phase the magnetizations can be expanded up to order (9(A~^/^) in 
terms of bosonic operators Oq,, aj^ : 

M^ = A/n-2aK, Ml = -%^¥l^{al-aa), M'^ = ^Nj^{ai+aa), (75) 
yielding the following expression, still bearing the Gallavotti-Cohen symmetry : 

^(A)= ^A(l-A)(/52 + a2)(5^(at +a„^^' 



2 

J]/5^((l - \)a\ + \a^){al + a„) + 2 ^ a^a^ + 2A/? 



n 

a, 7 



(76) 



where we define and recall that : 

^^ = -E^- P=-Y.P^ (77) 



n — ' n 

a a 

As in section in we find yu(A) is given by : 



1 /"dcj det(r + r]) 2 2x 



where F, F and Vl are 2n x 2n matrices defined by blocks in the following way : 

r(A) = 2 I ; - 1, F = 2 I ; - \,n=\ ^"^^ I (79) 





-iuld 



with elements : 



Aa, = -A(l - A) (/32 + ^2) - - {(3a + /3^) (80) 
n n 

= -A(l - A) {(3^ + ^2) - 1^ {(3a + (3,) (81) 
' n n 

Ca^ = -A(l - A) {(3^ + a^)-- ((1 - X)(3a + A/3-,) + 5^-^ (82) 
n n 
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Tedious yet straigtforward calculations lead to the following expressions : 

det(r + n) = (-l)"(cu2 + 4)"-2 {io^ + iol) {io^ + iol) (83) 

det(f + n) = (-1)"(C^' + 4)"-2 (^2 ^ ~2 ^) ^^2 ^ ~2 ^) (g4^) 

a;± = V(2 - /5)2 + 4A(1 - X)a^ ± ^(3^ + 4A(1 - A)^^ (85) 

(Jj± =2-(3 + 2X{l-X){l3'^ + a^)± 

^4A(1 - A)a2 + (2A(1 - A)(/?2 + ^2) - /5)2 (86) 

The expression of /i(A) simplifies into: 

^(A) = 2-f3- v/(2-/?)2 + 4A(l- A)a2 (87) 
We again obtain the full spectrum of the evolution operator: 

Sp HiX) = {2k + Lj4X)i + a7_(A)£' - /i(A)}M,^.eN (88) 

whose degeneracy is C^^^ • In the high temperature phase, the modes of the evolution 
operator split into two groups: two modes are similar to the n = 2 bath case, and 
n — 2 modes relax with constant rate 2. The slowest relaxation time of the master 
equation evolution operator is the same as that found for n = 2. The components of the 
probability current along the corresponding n — 2 directions cancel in the NESS, thus 
leaving the probability flow lines as ellipses, as is the case for n = 2. 

Finally, we can remark that the result ()87p for n{X) could also have been found 
by studying the set of Langevin equations associated to the Fokker-Planck equation of 
the system. However, the Langevin formalism does not render the Gallavotti-Cohen 
symmetry explicit at intermediate steps of the calculations, as opposed to the operator 
approach chosen in Sec. 4. 
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